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Abstract — This paper presents an algorithm for computing 
inner estimates of the regions of attraction of limit cycles of a 
nonlinear hybrid system. The basic procedure is: (1) compute the 
dynamics of the system transverse to the limit cycle; (2) from the 
linearization of the transverse dynamics construct a quadratic 
candidate Lyapunov function; (3) search for a new Lyapunov 
function verifying maximal regions of orbital stability via iterated 
of sum-of-squares programs. The construction of the transverse 
dynamics is novel, and valid for a broad class of nonlinear hybrid 
systems. The problem of stabilization of unstable limit cycles will 
also be addressed, and a solution given based on stabilization of 
the transverse linearization. 

I. Introduction 

Nonlinear dynamical systems exhibiting oscillating solu- 
tions are found in an extraordinary variety of engineering and 
scientific problems. For example: cellular signalling, radio- 
frequency circuits, walking robots, and population dynamics. 
Stability analysis of such oscillations has a long history, 
with the modern theory going back to Poincare. For planar 
systems substantial qualitative insight can be gained, however 
for higher-order systems almost all stability results are local. 
In this paper, we aim to characterize regions of stability to 
limit cycles of nonlinear systems. The problem of orbital 
stabilization via feedback conttol will also be addressed. 

A major motivation for the work in this paper is control of 
underactuated "dynamic walking" robots (|1|). These robots 
can exhibit efficient, naturalistic, and highly dynamic gaits. 
However, control design and stability analysis for such robots 
is a challenging task since their dynamics are intrinsically 
hybrid and highly nonlinear. Local stabilizing control design 
has been investigated via hybrid zero dynamics (0, 0) and 
transverse linearization (see, e.g., 0, 0, |[6lD. 

As well as being of interest in their own right, estimates of 
regions of stability will be an enabling technology for planning 
transitions among a library of stabilized walking gaits (Q), 
and for constructive control design and motion-planning ([8]). 

The most well-known tool for analysis of limit cycles is the 
Poincare map: orbital stability is characterized by stability of 
an associated "first-return map", describing the repeated passes 
of the system through a single transversal hypersurface. Often 
a linearization of the first-return map is computed numerically, 
and its eigenvalues can be used to verify local orbital stability. 
Since the system's evolution is only analyzed on a single 
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surface, regions of stability in the full state-space are difficult 
to evaluate via the Poincare map. 

A related technique known variably as "transverse coordi- 
nates" or "moving Poincare sections" also has a long history 
and was certainly known to exist by Poincare, however has not 
been much used in applications until recently due to difficulty 
in the relevant computations (see [9]). With this technique, a 
new coordinate system is defined on a family of transversal 
hypersurfaces which move about the orbit under study. In most 
cases, it is also used to study local stability, however as we 
will show it can be adapted to characterize regions of stability 
in the full state space. 

The construction of the transverse dynamics is novel, an- 
alytical, valid for a large class of systems, and amenable 
to analysis via sum-of-squares (SoS) programming. Previous 
work by the author and colleagues has utilized a construc- 
tion specifically for Lagrangian mechanical systems (see 0, 
0, 0). The new construction is also useful for design of 
stabilizing controllers for non-periodic trajectories of highly 
nonlinear systems without this special structure, e.g. in ifTUl 
a quadruped robot bounding over rough terrain was stabilized 
using a preliminary version of the construction in the present 
paper. 

In this paper we present an algorithm to compute a conser- 
vative estimate of the region of stability to a known periodic 
solution of hybrid nonlinear system. The method we propose 
is to construct the transverse dynamics in regions of the orbit, 
and then utilize the well-known SoS relaxation of polynomial 
positivity which is amenable to efficient computation via 
semidefinite programming (see, e.g., ifTTl . Ifl2l . Ifl3l ). The 
sum-of-squares relaxation has been previously used to char- 
acterize regions of stability of equilibria of nonlinear systems 
(see, e.g., HI, E2, Ml 

The method we propose has been tested on a number of ex- 
ample systems, presented in the companion paper IfTTl . There 
is comparatively little work on computing regions of stability 
of limit cycles. The proposed method has aspects in common 
with the surface Lyapunov functions proposed in 1 18 1, however 
that method was restricted to piecewise linear systems. The 
technique of cell-to-cell mapping, proposed by [19|, improves 
the efficiency of exhaustive grid-based methods of regional 
analysis and has been used in analysis of walking robots 
(1 20 ]), however the computational cost is still exponential in 
the dimension of the system. 
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II. Preliminaries and Problem Statements 

A. Stability of Limit Cycles 

Consider an autonomous dynamical system, which may be 
continuous or hybrid, with a state x £ R n . The solution, or 
flow of the system is denoted by $>(xq, t), i.e. x(t) = $(xo,t) 
is the solution at time t > of the dynamical system in 
question from an initial state x(0) = xq. If this solution exists 
and is unique then &(xo,t) is well-defined. Note that since 
the system is autonomous, we can consider flows starting from 
t = without any loss of generality. 

Suppose the system has a non-trivial T-periodic orbit, i.e. 
T > is the minimal period such that x* (t) = x* (t + T) for 
all t, and one would like to analyze the stability of this orbit. 
It is well-known that such a solution cannot be asymptotically 
stable in the standard sense, since perturbations in phase are 
persistent. The more appropriate notion is orbital stability. The 
definitions in this section are all standard (see, e.g., [9|, [2 1 J). 

Definition 1: Consider non-trivial T-periodic solution x* (t) 
of a dynamical system with flow $(•, •), and let T* denote 
the solution curve: T* = {x £ W n : 3t £ [0,T) : 
x = x*(t)}. The solution £*(•) is said to be asymptotically 
orbitally stable if there exists a b > such that for any Xo 
satisfying dist(xo, T*) < b the solution exists, is unique, and 
dist($(x ,f),r*) — !• as t -> oo. 

The distance to a set is defined in the usual way: dist(x, V*) = 
inf j, e r* \y — x\ with | • | the Euclidean norm in R n . 
A stronger statement is exponential orbital stability: 
Definition 2: A T-periodic solution x*(-) is said to be 
exponentially orbitally stable if it is orbitally stable and 
furthermore there exists a&>0,AT>0,c>0 such that 
for any xq satisfying dist(xo, T*) < b we have 

dist($(a;o,i)> r *) < K dist(x , r*)e _c *. 

The primary aim of this paper is to characterize regions of 
orbital stability: 

Definition 3: A set R C K." with non-empty interior R % and 
r* C R l is said to be an inner estimate of the region of stabil- 
ity of x*{-) if for all x £ R we have dist($(x ,t),r*) -> 
as t — > oo. 

B. Problem Statements 

This paper will suggest algorithms for three analysis and 
control problems for limit cycles: 

Problem 1: Given an autonomous smooth nonlinear system 
with state iel": 

x = f(x) (1) 

and a non-trivial T-periodic solution x* (t) = x*(t + T) of ([T| 
such that f(x*(t)) ^ for any t £ [0,T]. Characterize the 
stability of £*(•) and if it is exponentially orbitally stable then 
compute an inner estimate of its region of stability. □ 

Many systems are best modelled by a combination of 
smooth nonlinear dynamics with occasional moments of in- 
stantaneous change in the state. This may model a discrete 
change in a switching controller, or a period of extremely fast 
change in the system state, e.g. an impact event in a physical 
system. 



The second problem is to estimate regions of stability for a 
class of hybrid nonlinear systems: 

Problem 2: Consider an autonomous hybrid nonlinear sys- 
tem with state x £ M. n and switching dynamics defined 



between hyperplanes: 

x = f{x), x£S~ (2) 

x+ = A (a;), x G S~ . (3) 

Suppose /(•) and A(-) are smooth and A : S~ —> S + where 

S~ = {x : c'_x = d^,g{x) > 0}, (4) 

S + = {x: c' + x = d + }, (5) 

c_,C-|_ £ M. n , and d-,d+ £ K. Suppose x*(-) is a non- 



trivial T-periodic solution that undergoes N impacts at times 
{ii, ijv} + kT for integer k. We will assume that 
the impacts are not "grazing", i.e. d_f(x*(ti)) ^ and 
d + f{x*{ti)) £ for all i. 

The problem statement is to characterize the stability of 
) and if it is exponentially orbitally stable compute an 
inner estimate of its region of stability. □ 

For simplicity of expression we will consider the problem 
with a single switching surface and a single set of continuous 
dynamics, however the extension to multiple switches and 
multiple continuous phases is trivial. It is also straightforward 
to have the dimension of the continuous system change 
between different phases. 

The main practical restriction in this class is that switch- 
ing surfaces are planar. This is quite a strong restriction, 
but it greatly simplifies proving orbital stability via planar 
transversal surfaces, since we can make the transversal surface 
line up with the switching surfaces before and after impact. 
Furthermore, it is true for some important and common models 
of walking robots such as the rimless wheel and the compass 
gait. For other systems it may be possible to contstruct a 
change of coordinates such that the impact map is planar, e.g. 
for a multi-link walking robot one can choose one generalized 
coordinate to be the distance from the swing foot to the ground 
plane. 

The third problem is one of feedback orbital stabilization: 
Problem 3: Consider a controlled, possibly hybrid, system 
with a control input u £ K m : 

x = f(x,u), x £ S~ (6) 
x + = AO), x£ S-. (7) 

Suppose this system has a T-periodic solution £*(•) that un- 
dergoes impacts per period at times {ti,t2, tjy} + kT for 
integer k, and is generated by a piecewise continuous control 
signal, u*(t) = u*(t + T). If this solution is not exponentially 
orbitally stable then, if possible, construct a state-feedback 
exponentially orbitally stabilizing controller and compute an 
inner estimate of its region of orbital stability. □ 

III. Regions of Stability for Continuous Systems 

In this section we propose a solution to Problem [T] 
The process we propose for finding regions of orbital 
stability is based on the construction of a smooth local change 
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Fig. 1. Left: transversal surfaces S(t) around the target orbit x* , with a 
particular solution x (f ) converging to the orbit x*. Right: a Lyapunov function 
defined on a transversal surface. 



of coordinates x — > (x±,t). At each point t G [0, T] we define 
a hyperplane S(t), with 5(0) = 5(T), which is transversal to 
the solution T*, i.e. i*(t) g S(t). 

Given a point x nearby x*(-), the scalar r G [0, T) 
represents which of these transversal sufaces 5(r) the current 
state x inhabits; the vector x± G is the "transversal" 

state representing the location of x within the hyperplane 5(r), 
with xj_ = implying that x — x*(r). This is visualised in 
Figure [T] 

We will show that in some region around the target orbit r*, 
the dynamics in the new coordinate system are well-defined 
and have the form: 



Xj_ = A(t)x± + h(x±, r), 
t = l + g(x±,r), 



(8) 
(9) 



with h = 0(|xj_| 2 ) and g = 0(\x±\) near x± — for all r. 
The fact that there is a differential equation for r corresponds 
to the following fact: if the state x(ii) G S(t%), then after 
a short time interval St it does not necessarily follow that 
x{h+8 t ) G S(n + S t ). 

The main reason for this decomposition is that if we can 
prove that x± — > then we have proven orbital stability of 
x*. 

Associated with the above nonlinear system is the transverse 
linearization: a first-order approximation of the x± dynamics 
expressed as a periodic linear system: 



x± = A(t)x± 



(10) 



for t G [0,T). 



It is known that the periodic solution x* of the nonlinear 
system x = f(x) is exponentially orbitally stable if and only 
if the periodic linear system ( |T0] > is asymptotically and hence 
exponentially stable (see, e.g., (9), GTl ). 

The process for computing regions of orbital stability is as 
follows: 

1) Select a family of transversal surfaces 5(r), and the 
associated transformation x — > (x_l,t). 

2) Compute the nonlinear dynamics in this new coordinate 
system as well as a periodic linear system representing 
the dynamics of x± close to the orbit: the transverse 
linearization. 

3) Construct a candidate quadratic Lyapunov function as- 
sociated with the transverse linearization via standard 
techniques from linear control theory. 

4) Using this result as an initial seed, iteratively solve 
a sequence of sum-of-squares programs to compute 
maximal regions in which a Lyapunov function can be 
found verifying both well-posedness of the change of 
coordinates and orbital stability for the true nonlinear 
dynamics. 

The details of each step are given in the following four 
subsections 

A. Selection of a Set of Transversal Surfaces 

Suppose we have a periodic orbit x* (t) = x* (t + T) with 
x*(t) ^ OVt G [0, T). At each point x*(t) of the target orbit 
we define a transversal surface 5(r) in the following way: 

S(T) = {yeR n :z(T)'(y-x*(T))=0} 

where z(t) : [0, T) — > M. n is a smooth periodic vector function 
to be chosen. We will enforce that z(t) has bounded derivative. 
In the literature on the use of transversal coordinates to prove 
local properties about periodic solutions, it is common to 
choose z(t) = f(x*(r)) ||9], l22l . That is, the transversal 
planes are orthogonal to the current motion of the system. 

However, orthogonal transversal planes are often a bad 
choice when performing analysis on larger regions around the 
orbit, and allowing some freedom in z(t) is highly beneficial. 
The reason is, roughly speaking, that singularities will occur 
in the change of coordinates x — > (x± , r) near sections of 
x* with large curvature. This will be made more precise in 
Section |VT] 

The primary requirement is that the resulting 5(r) are still 
transversal to the orbit, which is guaranteed if there is some 
S > such that z(t)'/(x*(t)) > 5 for all r G [0,T). I.e., 
z(t) is never orthogonal to /(x*(r)). 

As a technical condition, we require that z(r) be Lipschitz 
on each continuous interval. For simplicity of derivations, 
and without loss of generality, we will further enforce that 
||z(r)|| = 1 for all r. If planes orthogonal to the system motion 
are desired, we can take z(t) = /(x*(t))/||/(:e*(t))||. 

B. Construction of a Moving Coordinate System 

Having chosen a set of hyperplanes, defined by z(t), we 
now construct a smoothly r-varying coordinate system upon 
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this subspace. For n = 2, construction of a basis is trivial: 
pick, e.g., [-Z2, zi]' as the basis vector for II. For n > 3, 
a constructive can be adapted from a method in E3l and (9] 
Ch. VI]. 

1) Choose a vector w such that tu is not collinear with z(t) 
for any r. 

2) Choose a fixed orthonormal basis rjj, j — 1,2, ...,n for 
K n with w as its first element. 

3) for each r € [0, T), define the plane containing both w 
and z(t), and the rotation matrix R(t) which takes w 
to be collinear with z(r) rotating in this plane. 

4) Define £,-(t) = R(T)r]j then £ 2 (V), £„(t) form a 
basis for the transversal coordinates. 



is 



(m + z(r)), j = 2,3,..,n. (11) 



Note that since r/i and z(r) are, by construction, unit vectors 
which are not collinear, ti' x z(t) < 1 for all r so ( fTT) is well- 
defined. 

The following lemma suggests that in practice the vector w 
in the above construction can be chosen at random: 

Lemma 1: Suppose z(t) : [0,T] — > M. n ,n > 3 is piecewise 
continuous with a finite number of points of discontinuity ij, 
and Lipschitz on each interval [ti,t i+1 ). Then the set of unit 
vectors w such that w = ±z(£) for some t £ [0,T] has 
measure zero. 

Proof: For each interval f,*+x) consider a closed interval 
[tj,£i_l_i] with z(ti + i) temporarily defined as 



z{U+i) 



t<ti 



lim 



z(t) 



Now, since z(f) is unit-length, we can consider the set z(t),t £ 
[ii, fj+i] to be a curve JE^ on the unit sphere R". Since has 
bounded first derivative, it is Lipschitz on [t$, and hence 
the curve Zj is rectifiable, and it is known that a rectifiable 
curve on a unit sphere in W 1 , n > 3 covers a set of finite 
measure. 

Since there is a finite number of continuous intervals, the 
union |J. Z{ is clearly also of zero measure. □ 

Having this basis defined, we also construct the projection 
operator: 



H(r) 



which defines the mapping x 
then 



Ur)'\ 

► (x±,t). That is, if x € S(t) 



x± = H(t)x. 

Note that a given x £ R" will in general be in more than one 
transversal plane, i.e. we can have x £ S(ti) and x £ S(t 2 ) 
with n 7^ T%, though this will not cause any problems for the 
proposed method. 



C. Transverse Dynamics and Linearization 

Theorem 1: The dynamics of the system in the new coor- 
dinates (x±,t) are given by: 

" d 

dr~ 

-n(T)/(z*(r))f, 
z(ryf(x*{T)+IL{T)'x ± ) 



Xj_ = 



: n(r) 



U(t)'x ± + n(r)f(x*(r) + n(r)'x ± ) 

(12) 
(13) 



z(r)'f(x*(r)) - *P- H(r)' 



x± 



Proof: Consider the transversal coordinate x±(t) = 
H(r)(y(t)—x*(T)) where r is such that z{r)'{y— x*{t)) = 0. 



An explicit formula for £j (t) in terms of T]±, rjj and z(t) x± = f 



dr 



n(r) 



[y(t)-x*(r)}+U(r) \f(y)-f(x*(r))t] 



Since y = x*(t) + n(r)'a;_L this can be written as ( fT2| 

To find the dynamics of r, consider the orthogonality 
condition: 



F(t,r) = z(T)'(x(t) - x*(t)) = 0. 



(14) 



Since this remains true under evolution of the system, the 



total derivative of F(t, r) is zero: F(t, r) 



dF , 



and, by the implicit function theorem, if ^ ^ we have 



dF 
di 



0. 



t = — (^p) %f-, after some straightforward manipulations 
from ( |14) , we find that 

z(r)'f(x(t)) 



z(T)'f(x*(T))-^(x(t)-x*(r))' 

As above, this can be given in terms of r and x± 

U(r)(x(t)-x*(r)): 

z(Tyf(x*(T)+U(r)'x ± ) 



(15) 



z(r)'/(^(r))-M^'n(r> 



Xj_ 



(16) 



□ 



7 J Transverse Linearization: The transverse linearization is 
the following (n — 1) -dimensional T-periodic linear system: 

i = A(i)z (17) 

where A(t) comes from differentiating ( fT2| ) with respect to 

x±: 

d 



A(t) = 



dt 



U(t) 



n(ty + u(t) df{x d ' x it)) n(ty 



-U(t)f(x*(t)) 



dt 



where 



dx± 



z{t) 



, df(x*(t)) 
dx 



my 



5f)'n(t)' 



z(t)'f(x*(t)) 



(18) 



(19) 



Remark 1: In two special cases some simplifications are 
possible: if the system is planar, i.e. n = 2, then 
[^n(r)]n(r)' = so ([12]) and (TH) can be simplified by 
removing the first term from each. If transversal planes orthog- 
onal to the system motion are chosen then II(t)/(x*(t)) = 
so ( p~2] > and ( fTS) can be simplified by removing the last term 
from each. 
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Remark 2: Note that in the selection of z(t) we have 
enforced that z(t)' f {x* (t)) > 5 for all r, so if x± remains 
sufficiently small we will have: 

«(t)7(jb*(t)) - U(tYx ± >a>0 

dr 

for some a, so the dynamics of r are well-defined. This 
condition breaks if there is a continuum of r such that x* (t) 
satisfies f(x*(T))'y — x*(t) = 0. E.g. for constant speed 
curves = 0) this means \y — x*(t)\ is exactly the radius of 
curvature of the the target orbit, i.e. l/\x\. In Section VI we 



will discuss optimizing z{t) so as to maximize the regions in 
which the dynamics of r are well-defined. 



D. Verification of Orbital Stability 

We now state conditions which guarantee regions of sta- 
bility, giving a solution to Problem [T] It will then be shown 
how to optimize regions for polynomial systems via a SoS 
relaxation. 

Theorem 2: Suppose there exists a function V : W 1 ^ 1 x 
K -> K such that D := {(x±,t) : V(x±,t) < 1} is compact 
and for which following inequalities hold for all (x±, r) 6 D: 



V(x ± ,t) > 0, x ± ^0(20) 
V(x±,t) < 0, x ± ^0(2l) 



dt 

V(0,t) = ~V(0,t) = 0, 

z(T)'f(x*(T)) - d ^ U(T)'X ± > 0. 
OT 



(22) 
(23) 



then the D is an inner estimate for the region of orbital stability 

Of £*(•). 

Proof: It follows from the smoothness assumptions on f(x) 
and ( |23~| that both x± and f are Lipschitz in D, and by the 
usual Lyapunov arguments ( |20) , pTj ) ensure that the compact 
set Z? is invariant. Therefore the solutions of r and from 
any initial conditions in D exist and are unique. Furthermore, 
( |2"T) , ( |22] i imply that from any initial condition in D we 
have x± — > 0, from which it follows that x* is orbitally 
asymptotically stable and D is an inner estimate of its region 
of stability. □ 

The optimization problem is then to search over Lyapunov 
functions V(x±,t) so as to maximize the size of the regions 
D such that these conditions can be verified. 

In practice, we will often strengthen the above statements 
by choosing some small constants Si > 0, 62 > 0, 63 > and 
requiring: 

Vix^^-S^xl 2 > 0, (24) 

^-V(x ± ,t) + S 2 \x\ 2 < 0, (25) 
dt 

z(Tyf(x*(r))- d ^U(T) , x x -5 3 > 0. (26) 

This improves the robustness of the numerical procedures for 
verification. 



1) Sums-of -Squares Relaxation: 

We now show how the regions can be computed using a 
SoS relaxation. For the purposes of this section, we assume 
that the vector field f(x) is polynomial in x. This being the 
case, let us fix a particular value of r and examine the formula 
for f: 



z(T)'f(x*(T)+IL(Tyx ± ) 



z(r)'f(x*(r)) 



dz{r) 



n{x±,r) 
d(x±, t) 



It is clear that both the numerator and denominator, n(x±,r) 
and d(x±,r) respectively, are both polynomial in x±. 

Furthermore, the well-posedness condition ensures that 
d(x± ,t) > 0, hence we can multiply both sides of condition 
( p5| ) by d(x±,r) without changing its validity, resulting in 
the condition DV(x±, r) < 0, where DV(x±,t) is given by 
( p7| ). It is straightforward to verify that this condition is also 
polynomial in x±. 

We can verify these conditions regionally using Lagrange 
multipliers l(x±) and m(x±) that are polynomial in x± with 
the following SoS constraints: 

- DV(x ± ,t) -1(x ± )(1-V(x ± ,t)) = SoS, (28) 

d(T,x 1 _)-5 2 -m{x 1 _){l-V{xt,T)) = SoS, (29) 

= SoS, (30) 

to(x_l) = SoS. (31) 

In practice we sample a sufficiently fine finite sequence 
Ti, i = 1,2, Ni such that t\ = 0, t^. = T. Then for each 
n and verify the conditions ([28) - pT| for each fixed r^. 

The objective is to maximize the regions satisfying ( |28] i - 
pi[ ). The decision parameters are the coefficients of V (within 
some restricted class) and the Lagrange multipliers I and m at 
each sample of r. This optimization is bilinear in the decision 
parameters, however a reasonable approach which has proven 
successful is to iterate between optimizing over multipliers and 
optimizing over V . 

If V is a polynomial of higher order than quadratic, then 
what exactly should be optimized in the search for V is open 
to some choice, but a natural candidate is to maximize the 
size of some ball in ||a;j_|| contained in the set V(x±, r) < 1. 
This can be expressed as a sum-of-squares program, and has 
been found to be effective for computing regions of stability 
for equilibria (see lfl4ll ). 



E. An Initial Candidate Lyapunov Function 

The iteration approach to solving a bilinear optimization 
problem requires a reasonably good initial guess. We pro- 
pose using the solution of the periodic Lyapunov differential 
equation for the transverse linearization as an initial guess for 
V(x x ,t). 

It is well known that the transverse linearization, a periodic 
linear system, is asymptotically (and hence exponentially) 
stable if and only if there exists a periodic quadratic Lyapunov 
function verifying stability, i.e. 

V{x±,t) = x'j_P(t)x ± 
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d d 
DV(x_l,t) := —V(x ± ,T)n{x ± ,T) + - — V(x±, r) 
OT ox± 



u(x±,t)t 



dr 



n(r) 



+d(x ± ,T)U(r)f(x*(T) + U{t)'x ± ) - U(r)f(x*(r))n(x ± ,T) 



n(r)'n 

+ d(.Tj_,T)c>i|:Z:jJ 2 . 



(27) 



with P(t) symmetric periodic positive definite (SPPD) such 
that 

P{t) + A{t)'P{t) + P(t)A(t) + H(t) < (32) 

for some SPPD H(t). One example is the unique periodic 
solution of the Lyapunov differential equation: 

P(t) + A(t)'P(t) + P(t)A{t) + Q(t) = (33) 

for any Q(t) = Q(t + T) > which can be freely chosen. 
In other cases, if ([TJ is the result of a linear control design 
procedure such as LQR or H°°, then P(t) may come from 
solution of an associated Riccati equation. 

Theorem 3: The following statements are equivalent 

1) x* is an exponentially orbitally stable T-periodic solu- 
tion of ([T}. 

2) The transverse linearization ( fT7j ) is exponentially stable. 

3) For each SPPD Q(t) there exists a unique SPPD so- 
lution P(t) to the periodic Lyapunov equation ( [33) . 
Furthermore, for a sufficiently small region around x* 
the Lyapunov function 

V(x) = 3> ± P(t)x u (34) 

verifies exponential orbital stability of x* . 
Proof: 

1^2: Consider dynamics in coordinates (x±,t). Due to 
exponential stability for a sufficiently small region around x* 
dynamics are well posed and Lipschitz. The linearization of 
the dynamics about the target orbit has the form 



d 


T 




'0 * 




T 


dt 


x±_ 




A(t)_ 




Xj_ 



(35) 



where * indicates a "don't care" element. This linearized 
system has a monodromy matrix (i.e. state transition matrix 
over the period T) of the form 



where ty± is the monodromy matrix of the transverse lin- 
earization ( fTTj i. It is known that if x* is exponentially orbitally 
stable, then has one eigenvalue Ai = 1 and all other 
eigenvalues satisfy |A,| < 1 (|9|). Due to the block diagonal 
form, this clearly implies that the eigenvalues of are 
all inside the unit circle, which implies the stability of the 
transverse linearization. 

2 «-» 3: Is a known result in periodic linear systems theory 
(see, e.g., d). 

3 — » 1: We will use Lyapunov's direct method with the 
Lyapunov function ( f34| ) to verify exponential stability. 



It is clear that V(x) = implies that x± = hence x € X* . 
Furthermore, from ([8]), d9), and ([32} we have 

V(x) < 2x±P{T)h{x±,T)+g(x x ,T)x' ± ~P{T) Xl _ 

dr 

—x±H(t)x± 

and h = 0(\x±\ 2 ) and g — 0(\x±\), so all terms except 
— x±H(t)x± are third-order in x±, therefore for sufficiently 
small we have V{x) < —a\\x±\\ for some a > 0. 

Now, since x* is a closed orbit, and by definition of x±, over 
some compact region around x* we have k\ dist(x, X*) < 
\\x±\\ < &2 dist(a;, X*) for some fci,fc2 > 0, hence 
&3 dist(x, X*) < V(x) < &4 dist(a;, X*) for some k%,ki > 
and V(x) < aki dist(a;, X*) which proves exponential orbital 
stability. □ 

A simple bisection search can be done to find a region 
V(x±,t) < p in which orbital stability is guaranteed. Then 
Vi(x±,t) — V(x±,t)/ p can be used as an initial seed in 
the search for more expressive Lyapunov functions around the 
orbit. 

IV. Hybrid Systems 

In this section, we extend the above results to the case of 
hybrid systems, and propose an algorithm for solving Problem 

m 

Suppose we find a vector function z(t) £ M™ for t 6 [0, T) 
such that 

z(t)'f(x*{t)) > 0,Vte[0,T) (36) 
z{t t ) = a iCi , V i = 1,2,. ..N, (37) 

for some sequence of scalars ai ^ 0. That is, z(t) aligns 
with the normals to the switching surfaces, and always makes 
a "small angle" with f(x*). The non-grazing condition on 
impact surfaces ensures that many such vector functions exist. 

With the coordinate system defined to line up with the 
switching surfaces, an entire region of the form V(x±, r) < 1 
for a time ti will map into the image of Si under A^. 

The switching update of the transversal coordinates is: 

x+ = II^A^T-r) +n(rn'xl) ~ x\t+)\ 
and the transverse linearization of the impact map is: 

A d = U(r+)^U(rry 

evaluated at x — x*(rf). 
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A. Stability Conditions 

For systems with impacts, stability conditions will be of the 
form 

jV(x±,t) < 0,(38) 

*(t)'/(z*(t)) - ^ HM'xx > 0,(39) 

(c'_ (** (rr ) + n(rf )'a:±) - d_) I. (xi , r) 

+ 5 (a;*(Tr) + n(r-)'x ± ) < 0,(40) 

for continuous phases over the regions V(x±,t) < 1 via 
Lagrange multipliers. Here I s (x±,t) is a Lagrange multiplier 
which is not constrained to be positive. The first constraint 
verifies stability via the decreasing Lyapunov function; the 
second verifies well-posedness of the change of variables; the 
third constraint ensures that the switching surface is not hit 
before it is expected. For some systems, it may be obvious 
that this will not happen, and the third condition could be 
dropped. At the switching times the conditions to be verified 
are: 

y(n(r+)[A l (^(rr)+n(rr)'^) 

-x*(t+)],t+) -V(x ± ,rr) < 0, (41) 
g(x*(Tr)+Il(Try x± ) > 0. (42) 

The first verifies stability, the second that the state is within the 
region of definition of the impact on the switching hyperplane. 
All constraints evaluated over the regions V(x±,t) < 1 via 
Lagrange multipliers. 

Theorem 4: Suppose conditions ([3~8"j)-(|4"2"]i hold in the re- 
gion D := {(x±,t) : V(x±,t) < 1} then D is an inner 
estimate of the region of attraction of the solution x*(-) of the 
hybrid system ([2]), ((3J. 

The proof is omitted because of space restrictions but is a 
minor variation on the proof of Theorem [2] and standard 
application of Lyapunov's direct method to impulsive systems 
(see, e.g., ED). 

1) Sum-of-Squares Relaxation: The above conditions for 
orbital stability of hybrid limit cycles can be relaxed to sum- 
of-squares programs analogously to (|2"8|i-(|3Tj) with additional 
SoS constraints added for the impact conditions ( pri) , |42|. 
Details are omitted due to space restrictions. 

B. Initial Candidate Lyapunov Function 

A natural candidate Lyapunov function for the hybrid case 
is the unique SPPD solution of the jump-Lyapunov function: 

-p = A'P + PA + Q, t^U 
P(rr) = A4WP(rf)MTi)+Qi, t = U 

A similar statement to Theorem [3] can be made applying the 
results of ll26l and (25). Details are omitted due to space 
restrictions. 



V. Controlled Systems and Orbital Stabilization 
Consider a controlled system: 

±(t) = /(&(*),«(*)) («) 

with x(t) G 1" and u(t) G R m . Suppose a nominal T- 
periodic solution has been found: x*(t), u*(t). We wish to 
study perturbations and orbital stabilization of this system. 
With the nominal control implemented as a function of time, 
the resulting perturbed system: 

x(t) = f(x*(t) + Ax(t),u(t)) 

is non-autonomous has extremely complicated dynamics and 
orbital stability cannot be characterized locally. 

In contrast, if u(t) is a feedback policy on x(t) then 
the region-of-attraction analysis reduces to the study of an 
autonomous system and the methods above can be applied 
directly. This analysis can of course be applied with any 
feedback control law, however the analysis of the transverse 
dynamics suggests a natural approach to the stabilization 
problem. 

Suppose a set of transversal surfaces is chosen, along with 
a projection f(-) : R n -> [0, T) such that x G S(f(x)) for all 
x in a neighbourhood of the solution £*(•), and we apply the 
control 

u(x) = u*(f(x)) + u. 
with u a stabilizing control, for example, 

u = -K{f{x))tt{f{x))[x - x*{f{x))] 

where K(-) : [0,T) -> R mxn is feedback gain. Note that the 
nominal command u* is being applied as feedback law based 
on the estimated phase t rather than as an open-loop function 
of time. 

This feedback gain K(-) can be designed using the con- 
trolled transverse linearization: 

± x =A(t)x±(t) + B(t)u(t), iG[0,T) (44) 

with A(t) is in ([18]) with f{x*{r) replaced by f(x*(r), u*(r)) 
and d ^ x d } T ^ replaced by d ^ x ^ u ^\ The formula for 
B(t) is 

B(r) = n(T) ^M^!M) _ n(r)/(/(r),^(r))^ 
where 

dr = 1 . , df( X *(T),U*(T)) 

du z(t)> f(x*(T),u*(T)) Z[T ' du 

Note that if orthogonal transversal surfaces are used then 
U(t) f (x* (t) , u* (t)) is zero and therefore the second term 
in ( |45] l is zero. 

Any K (•) which stabilizes the periodic system ( pB] ) locally 
exponentially stabilizes the target orbit. An associated Lya- 
punov function for the linear system can be used for regional 
verification of the closed-loop orbital stability as above. 

'This neighbourhood will be characterized by the computed regions of 
attraction of the closed-loop system, but is guaranteed to exist locally which 
is sufficient for the present arguments. 
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A. LQR via Jump-Riccati Equation 

As an example of possible control design procedure, we 
give the solution to the transverse-LQR problem via a jump 
Riccati equation. This solution is far from the only one, but 
is quite generic in that if the orbit is stabilizable in a certain 
sense, this method will always stabilize it. 

Theorem 5: The target orbit x*(-) of the original nonlinear 
system is exponentially orbitally stabilizable by smooth state 
feedback if and only if the associated transverse linearization is 
stabilizable. Furthermore, if this is the case then there exists a 
unique SPPD solution of the following jump-Riccati equation: 

-P = A'P + SA- PBR~ X B'P + Q, t^U 
P(rr) = A d fa)'P(rt)Mn) + Qi,t = t i 

with R, Q, Qi > 0. An associated locally stabilizing feedback 
controller for the original nonlinear system is then given by 

u(t,x±) = u*{t) - R- x B(t)'P{t)x±. 

Proof: The full proof is omitted due to space restrictions, 
however it follows the equivalence to exponential stabilization 
via smooth state feedback and stabilization of the transverse 
linearization from straightforward application of Theorem |4] 
The existence of the SPPD solution of the jump-Riccati 
equation follows from minor modification of the stabilizability 
properties of periodic impulsive linear systems ll26ll . Il24l . Note 
that the restriction that Q, Qi > can be relaxed somewhat 
to a discrete-continuous observability property of the cost 
function. □ 

VI. Optimization of Transversal Surfaces 

The method described above allows great flexibility in the 
choice of transversal surfaces, parameterized by their normal 
vectors z(t). In some cases, problem specific information 
might suggest a natural candidate (see examples). In others, 
the classical choice z(t) — /(x*(t))/||/(:e*(t))|| maybe suf- 
ficient. However, it is likely that for many practical examples 
some optimization of z(r) is appropriate. 

The well-posedness conditions break when 

z(r)'f(x*(r)) - TL<j)'x x = 0. 

It is natural to try to choose z(r) such that this is prevented 
from happening in a large region around the orbit, i.e. for a 
large ball in ||xj_||. 

Finding the smallest ||xj_|| for which this can happen is 
simply a linear-constrained least-squares problem, solved with 
the augmented cost function 

J(x_l,\) = y ±X± + \(z{T)'f(x*(r)) - ^ H(r) >x ± ) 
The usual computation leads to the following minimizer: 

t _ z(r)7(x*W)n(r)^ 

X ± 2 



and since we only care about the size of xj_, this simplifies to 





W(x*(t))| 




n(r)^ 





Note that ^V" * s orthogonal to z(t) and hence entirely in the 
row-space of LT(t), so |LT(t)^| = |^| and 

, t M , _ l*W(**(r))| 

I dr I 

making it unnecessary to compute II (r) at every step. 

The size of x ± (t) should be maximized in some sense. 
However, notice that the denominator has 9z S T ' as a factor. 
The optimum may have z{t) constant for some intervals of r 
which means |xj_(r)| would go to infinity, which may make an 
optimization difficult. It is therefore better-posed to minimize 
its inverse. E.g. one can minimize 

/ r r \ 1/p 

z{tY = argmin U wv^ j (46) 

s.t. z(r)7(x*(r)) > 5 (47) 
||z(r)| = 1 Vr e [0,T]. (48) 

for some small positive 8. The author has had success choosing 
p to be quite large, say between 10 and 100, so that it 
approximates an L°° norm but remains a smooth optimization 
problem. 

For systems with impacts, the integrand in the above opti- 
mization can be multiplied by a smooth shaping function </>(r) 
with 4>{Ti) = 0, since the transversal surfaces must be aligned 
with the switching surfaces and cannot be optimized. 

This is a nonconvex optimization, so it requires a good 
initial seed. A natural candidate optimization is z(t) = 
/(x*(t))/||/(x*(t))||. This initial guess always satisfies the 
trans vers ality conditions, but may not always satisfy the con- 
dition of alignment to switching surfaces. 

VII. Illustrative Examples 

The proposed method has been used by the author and 
colleagues to compute regions of stability for a number of 
example systems. The systems so far considered include: the 
van der Pol oscillator, which illustrates the importance of 
selection of z(t); the rimless wheel, a very simple model of 
an underactuated walking robot that illustrates the extension 
to hybrid systems; the compass-gait walker, a more complex 
model of underactuated walking with regions of stability that 
cannot be found in closed form, and illustrate the proposed 
method's use for nontrivial systems. A detailed discussion of 
these examples can be found in the companion paper ifTTl . 

VIII. Conclusions and Future Work 

We have presented a constructive procedure for computation 
of an inner (conservative) estimate of the region of attraction of 
a nonlinear hybrid limit cycle. Such a procedure was motivated 
by problems of control and motion planning of walking robots, 
but undoubtedly will find applications in many other fields 
where oscillating nonlinear systems are of interest. 
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In the case where a feasible but unstable periodic solution 
is known, the method of transverse coordinates can also be 
used to construct a feedback controller which guarantees 
exponential orbital stability. The regions of attraction of the 
closed-loop system can then be analyzed using the proposed 
technique. 

The method of analysis via Lyapunov functions and sum-of- 
squares relaxation has a natural extension to uncertain systems 
via storage functions and integral quadratic constraints (see, 
e.g., Il27l . Il28l ) or parameter-dependent Lyapunov function 
( 1 29 1). Robustness of the existence of periodic trajectories has 
previously been studied using an extension of the small-gain 
theorem in BUI . PP . Future work will include connections 
between such robust analysis tools and the present work on 
regions of stability. 
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